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As the realization of a fully operational quantum computer remains distant, quantum simulation, whereby one 
quantum system is engineered to simulate another, becomes a key goal of great practical importance. Here we 
report on a variational method exploiting the natural physics of cavity QED architectures to simulate strongly 
interacting quantum fields. Our scheme is broadly applicable to any architecture involving tunable and strongly 
nonlinear interactions with light; as an example, we demonstrate that existing cavity devices could simulate 
models of strongly interacting bosons. The scheme can be extended to simulate systems of entangled multicom- 
ponent fields, beyond the reach of existing classical simulation methods. 



Modelling interacting many-particle systems classically is 
a challenging yet tractable problem. However, in the quan- 
tum regime, it becomes rapidly intractable, owing to the dra- 
matic increase in the number of variables required to describe 
the system. Feynman [1] realized that an alternate approach 
would be to exploit quantum mechanics to carry out simu- 
lations beyond the reach of classical computers. This idea 
was the basis of Lloyd's simulation algorithm [2], a proce- 
dure for a digital quantum computer to simulate the dynamics 
of a strongly interacting quantum system. In contrast, there is 
also an analogue approach to quantum simulation, where the 
simulator's Hamiltonian is tailored to match that of the simu- 
lated system [3]. The complementary aspects of the analogue 
and digital methods, reviewed in [3-7], have led to a host of 
recent experiments [8-16]. 

To date, most experimental implementations of quantum 
simulation algorithms have been focus sed on the task of sim- 
ulating quantum lattice systems, with comparatively less at- 
tention paid to systems with continuous degrees of freedom. 
The archetypal example of a quantum system with a continu- 
ous degree of freedom is the quantum field. Currently, quan- 
tum simulations of quantum field theories have relied on dis- 
cretization of the dynamical degrees of freedom. One body 
of recent theoretical work is focussed on the analogue simula- 
tion of discretized quantum fields, using cold atoms in optical 
lattices [17-21] and coupled cavity arrays [22-24]. Comple- 
menting this are proposals for digital quantum simulation on 
a universal quantum computer of the zero-temperature [25] 
and thermal [26] dynamics of non-abelian gauge theories and, 
more recently, a digital quantum simulation [27, 28] of scat- 
tering processes of a discretized A(p 4 quantum field. 

In this paper we report on an analogue algorithm to simu- 
late the ground-state physics of a one-dimensional strongly 
interacting quantum field using the continuous output of a 
cavity-QED apparatus [29-33]. Our method involves no dis- 
cretization of the dynamical degrees of freedom; the simula- 
tion register is the continuous electromagnetic output mode 
of the cavity. The variational wave function generated in 
this way therefore belongs to an extremely expressive class, 



namely the class of continuous matrix product states, as we 
will show. We argue that our approach is already realizable 
with state-of-the-art cavity-QED technology. 

We concentrate on simulating quantum fields modelling 
collections of strongly interacting bosons in one spatial di- 
mension. These systems are compactly described in second 
quantization using the quantum field annihilation and creation 
operators ifi(x) and \ft(x), which obey the canonical commu- 
tation relations [ifr(x), &\y)] = S(x - y). The task is to deter- 
mine the ground state of a given field-theoretic Hamiltonian 
< K(0 r , (ft). The prototypical form of such a Hamiltonian is 

ik = J (f + W + N)dx (1) 

where f = d -^ d -M,W = J w(x - y)$Hx)p(y)$(y)$(x)dy, 
and N = -jiif/\x)ifr(x) describe the kinetic energy, two- 
particle interactions with potential w(x - y), and the chemi- 
cal potential, respectively. Our approach provides a quantum 
variational algorithm for finding the ground states of an arbi- 
trary Hamiltonian that is translation-invariant and consists of 
finite sums of polynomials of creation/annihilation operators 
and their derivatives. 

The apparatus proposed to simulate the ground-state 
physics of is a single-mode cavity coupled to the quantum 
degrees of freedom of some intracavity medium (Fig. 1); our 
proposal is not tied to the specific nature of the medium, so 
long as one or more tunable nonlinear interactions are present 
that are sufficiently strong at the single-photon level. Below 
we consider the example of a single trapped atom coupled to 
the cavity via electronic transitions. The system is described 
by a Hamiltonian H sys (A) that depends on a set of controllable 
parameters A, for example, externally applied fields. When 
the cavity is driven, either directly through one of its mirrors 
or indirectly through the medium, the intracavity field relaxes 
to a stationary state, and the cavity emits a steady-state beam 
of photons in a well-defined mode. 

The crucial idea underlying our proposal is to regard the 
steady- state cavity output as a continuous register recording a 
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FIG. 1. The output field E + (t) of a cavity-QED system is identified with a bosonic quantum field i(/(x). Since optical detection schemes 
correspond to expectation values of quantum-field operators, can be estimated via independent measurements of the cavity field. For 
example, the operators N, W, and T of Eq. (1) are determined, respectively, from measurements of (a) the output-field intensity, (b) Hanbury 
Brown and Twiss correlations, and (c) an interferometer with variable path length. 



variational quantum state ^(A)) of a one-dimensional quan- 
tum field with control parameters A as the variational parame- 
ters. This representation is chosen so that the spatial location 
x of the simulated translation-invariant field is identified with 
the value of the time- stationary cavity output mode exiting the 
cavity at time t = x/s. The arbitrary scaling parameter s is in- 
cluded in the set of variational parameters A. We complete 
this identification by equating the annihilation operator if/(x) 
of the simulated quantum field with the field operator E + (t) 
for the positive-frequency electric field of the cavity output 
mode [34], via fi(x) = E + (t)/ a/s. 

Recall that the variational method proceeds by minimiz- 
ing the average energy density of the variational state f(A) = 
0F(/l)|f + W + #|*F(/l)) over the variational parameters A. A 
key point in our scheme is that — with the identification of 
the field operators E + (t) and iff(x) in hand — the value of f(A) 
can be determined from standard optical measurements on the 
cavity output field, namely the measurement of Glauber cor- 
relation functions [35, 36], see Fig. 1. This result is easily 
seen for the Hamiltonian of Eq. (1). Thanks to the linearity 
of the expectation value, we can separately measure (f), (W), 
and (N). The expectation value of the chemical potential term 
corresponds to a function of the intensity of the output beam 
via (N) = -- s {E~(t)E + (t)). The kinetic energy term (f) cor- 
responds to the limit 

<f>= lim -^(g(X>(f + e l9 t + e z )'g (1 \t + e l9 t) 

ei,e2->0 S 5 E\€2 V 

- g V\t,t + e 2 ) + gV(t,t)), 

where g^\ti,t 2 ) - (E~ '{t\)E* (t 2 )); this quantity can be es- 
timated by choosing a finite but small value for e\ and e 2 . 
Note that this procedure does not amount to a simple space 
discretization because the output is a continuous quantum 
register. The final term (W) depends on two-point spatial 
correlation functions ffi(x)ffi(y){ff(y){ff(x)), which translate 



to measurements of g (2) (t 1 ,t 2 ) = (E~(ti)E~(t 2 )E + (t 2 )E + (ti)). 
The detection schemes to estimate all terms in the show- 
case Hamiltonian (1) are presented in Fig. 1. From a 
wider perspective, any Glauber correlation function g^ n,m ^ = 
(E~(h) • • • E-(t n )E + (t' m ) • • • E + (t[)), i.e., any n + m-point field 
correlation function composed of n creation operators E~ and 
m annihilation operators E + [35, 36], can be measured with 
similar, albeit more complex setups, such as in [37, 38]. Thus, 
upon identifying E~(t) and E + (t) with ifs\x) and ^(x), respec- 
tively, our scheme admits the measurement of any equiva- 
lent energy density oc (^(xi) • • • i(r\x n )i(r(x' m ) • • • tfr(x[)), and 
therefore ultimately the simulation of arbitrary Hamiltonians 

Once f(A) = < X F(^)| < K| X F(^)) has been experimentally es- 
timated for a given A, the next step is to apply the varia- 
tional method to minimize f(A). Minimization is carried out 
by adaptively tuning the parameters A in the system Hamilto- 
nian H sys (A) and iteratively reducing f(A), for example, using 
a standard numerical gradient-descent method. Once the op- 
timum choice of A is found, the resulting cavity output field is 
a variational approximation to the ground state of *H, and rel- 
evant observables of the field theory can be directly measured 
using the detection schemes of Fig. 1. We emphasize that our 
method applies also to cases where a numerical estimation of 
f(A) cannot be performed efficiently due to the size and com- 
plexity of the system to be simulated, and we suggest that this 
is exactly the strength of our approach. Moreover, the opti- 
mization may be performed experimentally without theoreti- 
cally calculating the cavity-QED system dynamics; indeed, it 
is not necessary to accurately characterize H sys or its relation 
to the adjustable parameters A. 

Why should the stationary output of a cavity-QED appara- 
tus be an expressive class capable of capturing the ground- 
state physics of strongly interacting fields? It is possi- 
ble to show that such states are of continuous matrix prod- 
uct state (cMPS) type, a variational class of quantum field 
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FIG. 2. Two-particle correlations in the Lieb-Liniger model are 
reproduced in simulations of an ion-trap cavity experiment. This 
critical model exhibits a transition between the superfluid regime 
for v « and the Tonks -Girardeau regime for v » 0, which is 
seen in the value of the correlation function at t = 0. (a) The 
Lieb-Liniger ground state is simulated for interaction strengths v = 
{0.07 (red), 3.95 (orange), 60.20 (yellow), 625.95 (green)}, and corre- 
lation functions ffi (0)ffi (x)fi(x)iff(0)) calculated as in [39, 40] us- 
ing 338 variational parameters, (b) Two-photon correlation functions 
g (2) (0 for an experiment with the parameters of [53]. Although there 
are visible differences, with just three variational parameters {g, £2, s] 
the transition in the correlation functions is approximately repro- 
duced. It is worth emphasizing how unusual it is for a variational 
calculation with only a few parameters to reproduce anything more 
than the coarsest features of a correlation function, e.g., if mean-field 
theory is used one does not obtain nontrivial correlators. Strikingly, 
the transition in (a) is captured even in the presence of realistic de- 
cay channels (inset). Note that this transition is analogous to that 
observed in [54]. 



states recently introduced for the classical simulation of both 
nonrelativistic and relativistic quantum fields [39-41, 65]. 
These states are a generalisation of matrix product states 
(MPS) [43-46], which have enjoyed unparalleled success in 
the study of strongly correlated phenomena in one dimension 
in conjunction with the density matrix renormalization group 
(DMRG) [47, 48]. It turns out that all quantum field states 
admit a cMPS description, providing a compelling argument 
for their utility as a variational class [49]. Crucially, the cMPS 
formalism turns out to be identical to the input-output formal- 
ism of cavity QED [50]. This identification was anticipated 
in [39, 40, 51, 52], and we elucidate it further in the supple- 
mental material. It implies that quantum field states emerging 
from a cavity are ofcMPS type and thus fulfill the necessary 
conditions for being a suitable and expressive class of varia- 
tional quantum states. 

Even though cavity-QED output states are of cMPS type, 
can a realistic system in the presence of decoherence repro- 
duce the relevant physics of a strongly interacting quantum 
field? As a test case, we demonstrate that the paradigmatic 
cavity-QED system, comprising a single trapped atom cou- 
pled to a high-finesse cavity mode, is capable of simulating the 
ground-state physics of an equally paradigmatic field, namely, 
the Lieb-Liniger model [55]. This model describes hard-core 
bosons with a delta-function interaction and is given by Eq. 
(1) with w(x-y) = vS(x-y), where v describes the interaction 
strength. Our simulator consists of a two-level atom inter- 



acting with one cavity mode, described (in a suitable rotating 
frame) by the on-resonance Jaynes-Cummings Hamiltonian 

H sys = g(a + d + + Q(<x + + &'), (2) 

where <x + is the atomic raising operator and a is the cavity 
photon annihilation operator, g the atom-cavity coupling, and 
Q the laser drive. The cavity-QED Hamiltonian H sys can be 
realized in various experimental architectures [29-33]. Here 
we choose the example of a trapped calcium ion in an optical 
cavity, with which tunable photon statistics have previously 
been demonstrated [54], and we show in the supplemental in- 
formation (see below) how g and Q can function as variational 
parameters. 

In an experiment, to measure the variational energy den- 
sity f(A), the output beam would be allowed to relax to steady 
state, the intensity /, and g (2) functions estimated as de- 
picted in Fig. 1, and f(A) finally estimated from postprocess- 
ing this data. Measurement schemes 2a and 2b of Fig. 1 
are just the standard laboratory techniques of photon detec- 
tion and Hanbury-Brown and Twiss interferometry. Measure- 
ment scheme 2c represents an interferometer with variable 
path length that is used to estimate the derivative of the quan- 
tum field in the kinetic energy term (f). Length shifts on the 
mm scale correspond to ps values of e\ and 62 in the estima- 
tion of (t); as these values are six orders of magnitude smaller 
than the relevant timescales of the experiment, they are con- 
sidered sufficiently small. 

Obviously the test model chosen here is simple enough to 
admit a classical simulation, which we carry out for the ex- 
perimental parameters of [53]. Exploiting a simple gradient- 
descent method, we find the values of A = {g, Q, s} that min- 
imize f(A) for a given value of v. This procedure is repeated 
over a range of v values of interest, and the corresponding 
optimized values of A are then used to compute quantities of 
interest, e.g., spatial-correlation functions as shown in Fig. 2. 
Remarkably, we find that just these three variational param- 
eters A = (g, Q, s), when varied in the experimentally feasi- 
ble parameter regime of [53, 54] in the presence of losses, 
allow for a quantum simulation of Lieb-Liniger ground-state 
physics. One expects that the ground-state approximation 
would improve by increasing the dimension of the auxiliary 
system and by allowing sufficiently general internal couplings 
and couplings to the field. In the context of atom-cavity sys- 
tems, this can be done by making use of the rich level structure 
of atoms (i.e., Zeeman splittings) and making use of motional 
degrees of freedom. For sufficiently complex intracavity dy- 
namics a classical simulation will become unfeasible, and at 
the same time it becomes conceivable that such a simulator 
will outperform the best classical methods. 

The reliability of a quantum simulator can be compromised 
by decoherence effects, as was recently emphasized in [56]. 
Our simulation of the ion-cavity system includes both cav- 
ity decay (at rate a:) and decay of the ion due to spontaneous 
emission (at rate y). The cavity decay rate rescales the pa- 
rameter s linking measurement time and simulated space, and 
thus it can be considered as a variational parameter itself. We 
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emphasize that cavity decay does not function as a decoher- 
ence channel in our scheme but is rather an essential element 
of the cMPS formalism. In contrast, spontaneous emission 
sets the limit for the timescale over which coherent dynamics 
can be observed. For our present example, we can show that 
the regime of strong cooperativity C = g 2 I icy > 1 is suffi- 
cient to allow a simulation of the Lieb-Liniger model despite 
detrimental losses. The Lieb-Liniger model exhibits nontriv- 
ial variations (Friedel oscillations) of the two-point correla- 
tion function on a length scale £ = {{//(x)^ (x))~ l . This length 
scale in the simulated model translates to a time scale over 
which the stationary output field (as the simulating register) 
should exhibit similar nontrivial features. From the previously 
established scaling rule one finds that the required time scale 
is r = g/s = (E~(t)E + (t))~ l . By means of the cavity input- 
output relations, we relate the output photon flux to the mean 
intracavity photon number (E~(t)E + (t)) = K(cfta). In the bad 
cavity limit k » g the cavity mode can be adiabatically elim- 
inated. In this case (a? a) ^ (g/tf) 2 , such that r = K/g 2 . On 
the other hand, the characteristic decoherence time of the ion 
is determined by the spontaneous emission rate 2y. Beyond 
a time l/2y the second order correlation function g (2) will be 
trivial. We therefore demand r < l/2y, which is equivalent 
to the requirement of a large cooperativity C > 1. For the 
exemplary case of the ion-cavity system considered above the 
cooperativity was indeed C - 1.8, see supplementary mate- 
rial. While equivalent conditions must be determined on a 
case-to-case basis, we expect that nontrivial quantum simula- 
tions in cavity QED will not be possible in the weak-coupling 
regime. Finally, there are overall losses associated with scat- 
tering and absorption in cavity mirrors, detection path optics, 
and photon-counter efficiency. However, while these losses 
reduce the efficiency with which photon correlations are de- 
tected, they do not otherwise affect the system dynamics. 

A natural question is when our scheme would provide a 
practical advantage over classical computers in the simulation 
of quantum fields. We expect this to be the case in particu- 
lar for the simulation of fields with multiple components, or 
species of particles, with canonical field-annihilation opera- 
tors a = 1,2, . . . , N. This situation arises in at least 
two settings: firstly, for vector bosons in gauge theories with 
gauge group SU(N), and secondly, in the nonrelativistic set- 
ting of cold atomic gases with multiple species. Variational 
calculations using cMPS fail in these settings, as the number 
of variational parameters must scale as D ~ 2 N . On the other 
hand, in a cavity-QED quantum simulation multiple output 
fields are naturally accessible via polarization or higher or- 
der cavity modes, and at the same time large effective Hilbert 
space dimensions can be achieved, e.g., with trapped ions or 
atoms. With N > 10, substantial practical speedups are al- 
ready expected with respect to the classical cMPS algorithm, 
which requires a number of operations scaling as 2 3xN . 

The realisation that ground- state cMPS and the field states 
emerging from a cavity are connected can be exploited to char- 
acterise the correlations of the emitted light. Indeed, we obtain 
a simple criteria to determine when the correlations in the light 



field are nonclassical: if it turns out that the simulated hamil- 
tonian is quadratic in the field operators; and (b) contains only 
"ultralocal" terms, i.e., no derivatives in the field operators, 
then the ground state is a trivial (i.e., gaussian) product state, 
and there would be no nonclassical correlations in the emitted 
light. 

The output of a cavity-QED apparatus admits a natural in- 
terpretation as a variational class of quantum-field states. We 
have demonstrated that this allows an analogue quantum sim- 
ulation procedure for strongly correlated physics using current 
technology. This result opens up an entirely new perspective 
for all cavity-QED systems which exhibit sufficiently strong 
nonlinearities at the single-photon level. This includes not 
only optical cavities coupled to atoms, but also superconduct- 
ing circuits with super-strong coupling to solid state quantum 
systems [57], as well as other nonlinear systems achieving a 
high optical depth without cavities, such as atomic ensembles 
exhibiting Rydberg blockade [58], or coupled to nanophotonic 
waveguides [59, 60]. Looking further afield, since the input- 
output and cMPS formalisms generalize in a natural way to 
fermionic settings [61-63], our simulation procedure might 
be applicable to cavity-like microelectronic settings involving 
fermionic degrees of freedom. We hope our work will inspire 
explorations of these promising directions. 
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Supplementary Material 
Continuous Matrix Product States and Cavity QED 

The matrix product state formalism has recently been gen- 
eralized to the setting of quantum fields in [39, 64, 65] giving 
rise to continuous matrix product states (cMPS). These states 
refer to one-dimensional bosonic fields with annihilation and 
creation operators $r a (x) and tfrl(x) (such that [^(x), t^(x)] = 
dapd(x - y)), and are defined as 

m = tro^lQ), (3) 
where |Q) denotes the vacuum state of the quantum field and 

U = V exp | J dxH cM ps (*) j • 

Here V denotes path ordering, and the (non Hermitian) Hamil- 
tonian is 

# cM ps (*) = G®1 + X^ ^ W W 

a 
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with Q and R a being D x D matrices acting on an auxiliary 
system of dimension D. tv aux in Eq. (3) is the trace over this 
auxiliary system. 

As was shown in [39, 40] an equivalent representation of 
the state |^) as defined in Eq. (3) is given by 

pcx lim twI^xo^OIQXQI^I^I^Uo,^)) (5) 

Xq,X\ — >oo I y 

where 

U(x , xi ) = 9 exp | J dxH cMFS (*) j , 

and |^) is an arbitrary state of the auxiliary system, 
which plays the role of a boundary condition at xo. The 
two states in (3) and (5) are equivalent in the sense 
that they give rise to identical expectation values for 
normal- and position-ordered expressions of field operators 
#(yi) . . . ft(y n )$(yn + i) • • • fom)), see [39, 40]. Note that the 
trace in (5) is a proper partial trace over the auxiliary system 
(in contrast to the trace in (3)). 

Consider, on the other hand, a cavity with several relevant 
modes described by annihilation/creation operators a a , 

(such that [a a ,cip] = 8 a p) which are coupled to some in- 
tracavity medium via a system Hamiltonian H sys . More- 
over, each cavity is coupled to a continuum of field modes 
([a a (a)),a^(a))] = S(oj - 6))5 a p) through one of its end mir- 
rors. (The generalization to the case of double- sided cavities, 
or ring cavities with several outputs per mode is immediate.) 
The total Hamiltonian for the cavities, the intracavity medium, 
and the outside field is 

H c QE D (t) = Hsys ® 1 

This Hamiltonian is written in an interaction picture with re- 
spect to the free energy of the continuous fields, and it is taken 
in a frame rotating at the resonance frequencies of the cav- 
ities. In the interaction picture and rotating frame each in- 
tegral extends over a band width of frequencies around the 
respective cavity frequencies. In the optical domain the Born- 
Markov approximation, which assumes K a {cS) = const, in the 
relevant band width, holds to an excellent degree. In this 
case it is common to define time-dependent operators E*(t) = 
J doj/ ' ^2jra a (a))exp(-io)t) 9 which fulfill [E+(t),E~(t')] = 
S(t - t'). As in the main text, these operators correspond to 
the electric field components at the cavity output, and are de- 
fined such that (E~(t)E + (t)) denotes the flux of photons per 
second. The Hamiltonian then is 

HcqedO) = H sys ® 1 + i^j (aa ® E~(t) - h.c.) . 

a 

If the outside field modes are in the vacuum, standard quan- 
tum optical calculations [66] show that this Hamiltonian is 



equivalent to an effective non-Hermitian Hamiltonian 

HcQEuit) = H S ys ® 1 - / ^ ^a\a a <g> 1 + i ^ <faa a E~(t). 

a a 

When this is compared to Eq. (4) the identification of the for- 
malism of cMPS and cavity QED is immediate. If we assume 
that the field and the cavity system is in a state |Q}®|0) at some 
initial time to the final state of the field outside the cavity at 
time t\ is given by expression (5), when we identify x = s t, 
and ^cMPs(^) = -iH cQ E D (x/s), that is 

Hx)=^-E + a (tl R a = A [^a a , Q = - l -H sys -V^ 

y a 

(6) 

with an arbitrary scaling factor s. Therefore, the state of the 
output modes of a cavity is always a continuous matrix prod- 
uct state. Formally these states are cMPS with an infinite- 
dimensional auxiliary system D — > oo. However, due to en- 
ergy constraints, the dimensions of the cavity system are ef- 
fectively finite. The relevant dimension of the cavity Hilbert 
space then sets the dimension D of the auxiliary system in the 
cMPS formalism. 

Quantum Variational Algorithm 

As an exemplary test case we demonstrated that the cavity 
QED system comprising a single trapped atom strongly cou- 
pled to a single high-finesse cavity mode is capable of simulat- 
ing the ground-state physics of the Lieb-Liniger model. The 
atom-cavity system is described (in a suitable rotating frame) 
by the on-resonance Jaynes-Cummings Hamiltonian 

H sys = g(d- + a + fr'tf) + D(o- + + <x~), (7) 

where & + is the atomic raising operator and a is the cavity- 
photon annihilation operator, g the atom cavity coupling, and 
Q the laser drive. Photons leak out of the cavity with leakage 
rate /c, and it is assumed that, in a real experiment, this output 
light can be measured by various detection setups. 
The Lieb-Liniger model, with Hamiltonian 

(f + W + N)dx 

oo 

J -co ctx (XX 

- J2^ (x)ijr(x)^dx, 

describes hard-core bosons with contact interaction of 
strength v. We performed variational optimizations for a 
range of values for v. We did this using a simple gradient- 
descent minimization of the average energy density f(A) = 
0F(/l)|f + W + N^iA)). Using the parameters that minimize 
f(A) we then calculate other quantities of interest, such as cor- 
relation functions for the simulated ground state field. We 
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have focussed on a gradient-descent algorithm for clarity; in 
practice a more sophisticated optimization procedure using 
the time-dependent variational principle, or conjugate gradi- 
ents, could be used. 

Our variational parameters A - (g, Q, s) enter f(A) as fol- 
lows. The expectation value of the energy density in the Lieb- 
Liniger model is 

f(A;v,p) = (f) + (W) + (N) 

where each term may be written in terms of the experimentally 
observed correlation functions via the correspondence $r(x) = 
E + (t)/ V?, giving 

(f > = lim (g«\t + 61, t + 6 2 ) - g (1) (* + € U t) 

-g (1) (M + e 2 ) + s (1) (U)), 
<#)= 4^0, 

<^> = -V ) a,o, 

where S^fo,^) = {E~ {h)E + {t 2 )) and g (2) (^ 2 ) = 
{E~(t\)E~(t2)E + (t2)E + (t\)). In an experimental simulation, 
these correlation functions would be measured directly in the 
laboratory, and the results would be fed back into a classical 
computer performing the optimization algorithm. However, 
for the purposes of our proof-of-principle simulation, we cal- 
culate the correlation functions directly by means of the input- 
output formalism. Making use of the cavity input-output re- 
lation E + (t) = EtJt) + ^a(t), where E^M) denotes the 
field impinging on the cavity at time t (which is assumed to 
be in the vacuum state), the expectation value may be written 
[39,40] 

f(A; v, Ai) = tr {([& R]f [Q, R] p ss } + v tr {(#) 2 R 2 p ss } 

-fitr{&Rp ss } 

where R and Q are defined in Eq. (6), and p ss is the unique 
steady state of the atom-cavity system, satisfying 

^ = -i [# sys , p ss ] + Kap^tf - ^ [a 1 ^, p ss ] = 0. 

Note that above, we write f(A;v,ju) to highlight the depen- 
dence of / on v and [i. Hereafter we set [i = 1 and minimize 
/(/i;v,yu) for a range of different values of v. Solutions for 
other values of /d can be obtained by means of a scaling trans- 
formation, as described in [39, 40]. 

The general outline of the algorithm to determine our opti- 
mum values of the variational parameters A, for a given choice 
of v and yu, is thus: 

initialize tol (tolerance), and e (step size) 
initialize A to an arbitrary value 
repeat 



set X = A 

calculate V/(A; v,yu) 

update experimental parameters as A <— A - eVf(X, v,yu) 
calculate f(X;v,jj) and f(A; v,yu) 
until \ f(A; v,/j) - f(X;v,/j)\ < tol 

Note that in our simulation, V/(/i;v,/i) is estimated numeri- 
cally by evaluating f(A + A ; ; v,//) for small values of A/, while 
in an experiment, Vf(X, v,ji) is found with the aid of measure- 
ments of <f), (W), and (N). 
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